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Abstract 


Science Applications International Corporation is currently developing the Electric 
Propulsion Interactions Code, EPIC, as part of a project sponsored by the Space Environments 
and Effects Program at the NASA Marshall Space Flight Center. Now in its second year of 
development, EPIC is an interactive computer tool that allows the construction of a 3-D 
spacecraft model, and the assessment of a variety of interactions between its subsystems and the 
plume from an electric thruster. These interactions may include erosion of surfaces due to 
sputtering and re-deposition of sputtered materials, surface heating, torque on the spacecraft, and 
changes in surface properties due to erosion and deposition. This paper describes the overall 
capability of EPIC and provides an outline of the physics and algorithms that comprise many of 
its computational modules. 


Introduction 


More than one hundred satellites are currently operating in space with onboard electric 
propulsion (EP). Many more missions with EP systems - both near-Earth and deep space - are 
being considered. A critical engineering issue in the integration of EP on spacecraft is the 
potentially unfavorable interaction of their high-energy plumes with surrounding critical 
components and diagnostic equipment. Such interactions may affect mission lifetime and 
sometimes even threaten mission success. As NASA considers more powerful EP technologies 
for long-duration missions, the need to predict pertinent interactions becomes even more critical. 
Only a few programs with the goal to produce large-scale, 3-D, global computer tools for the 
assessment of such interactions exist - two in the U.S. 1 ’ 2 and one in Europe. 3 One of the two U.S. 
programs aims in the development EPIC. 

EPIC is currently developed by Science Applications International Corporation (SAIC) as 
part of a 2-year effort that is sponsored by the Space Environment and Effects (SEE) Program at 


* Major portions of this paper have been presented at the AIAA 38 th and 39 th Joint Propulsion Conferences (AIAA 
paper numbers 2002-3667 and 2003-4871) 



the NASA Marshall Space Flight Center (MSFC). The objective is to produce an integrated 
modeling package that can be used to aid the design of spacecraft with onboard electric 
propulsion. The suite of computer tools is intended to allow both expert and novice users to 
investigate EP plume-S/C interactions. The code is built on the Module Integrator and Rule- 
based Intelligent Analytic Database (MIRIAD) architecture, which has formed the core of 
several NASA and Department of Defense programs. 

Overview 


In general, the spacecraft designer provides the following input to EPIC : satellite geometry 
and surface materials, thruster locations and plume parameters, case study parameters such as 
sputter yield coefficients, orbit and hours of thruster operation. The output may be in the form of 
contour plots of the plume map in space and of surface interactions on the 3-D spacecraft, 1-D 
plots along surfaces (e.g. erosion depth on a solar array as a function of distance from the 
thruster), integrated results over duration of mission (e.g. total induced torque in a given 
direction, total deposition of eroded material at a specific location on the spacecraft), list of 
results in text format for post-processing. EPIC ' s various components have already been used to 
asses EP-spacecraft interactions for both the government 1,4 and industry. 5 

The computer tool (built for the Win32 platform) consists of the following main components: 

• An enhanced version of the 3-D geometry-definition tool Object Toolkit (OTk) to 
facilitate interactive generation of the spacecraft geometry and materials. 

• A 3-D interactions code built using the Module Integrator and Rule-based Intelligent 
Analytic Database (MIRIAD) architecture. 

• A 2-D plume code that generates the plume map from an electrostatic thruster in the R-Z 
plane. A fast, 1-D Hall effect thruster (HET) algorithm is also part of the plume modeling 
capability. The algorithm is intended to provide a more quantitative estimate of the 
conditions near the thruster exit. 

EPIC_ Architecture 

The bulk of EPIC ' s capability is the 3-D plume-surface interactions code. The code is built 
on the MIRIAD architecture, which has formed the core of several NASA and DoD programs, 
and has allowed spacecraft designers to quickly and affordably model interactions between the 
spacecraft and the space environment. MIRIAD’s open architecture provides a framework for 
integrating a variety of physical models and their constituent data into a single executable 
application. The application then allows the user to define systems of interest and perform 
parametric studies with relative ease. Ideally, any variable in MIRIAD can be included in a 
parametric study. 

The main EPIC user interface and the code that integrates the various components are 
programmed in C# (C-sharp). C# is a new language developed by Microsoft that provides an 
easy-to-use object orientated way to program the windows user interface, and makes it extremely 
easy to combine components written in different languages using either the new Common 



Language Runtime (CLR) or the Component Object Model (COM) that has been the backbone 
of windows programs for years. The combination of C# and the CLR foundation classes is a very 
powerful and easy way to get the desired functionality in the user interface. 


The most natural and powerful way to communicate with the MIRIAD data server that 
coordinates EPIC calculations is through COM. It was therefore particularly important to be able 
to use COM components in EPIC. Any COM or dot net client can access EPIC ' s COM data 
server. For example EXCEL can be used to script calculations and manipulate data produced by 
EPIC. C# allows easy use of our COM- based 3-D viewer (MxOpenGlArena) for viewing the 
system and displaying surface information and spatial data. 

Main Computation Components of EPIC 


Spacecraft model 

Object Toolkit is EPIC ' s 3-D object generator. Originally developed to create spacecraft 
surface models for the Nascap-2k plasma interactions code, 6 OTk has been generalized to allow 
for the definition of spacecraft models for other analysis codes. In particular for EPIC, 
employment of “Special Objects” has been allowed to facilitate the specification of thrusters as 
illustrated in Figure 1. “Special Objects” are objects understood by the target application that are 
included in Object Toolkit primarily to specify and display their location. They are therefore not 
a part of the spacecraft surface model. Each special object instance has an object type name ( e.g . 
“Hall Thruster”), an instance name (e.g. “SPT100-1”), and a set of properties that includes, at 
minimum, its position (“x”, “y”, “z”), thrust unit vector, and the color with which it is to be 
displayed (“r”, “g”, “b”). The properties can be edited for each special object instance in similar 
fashion to the attribute property values. Figure 2 shows the OTk spacecraft model of Figure 1 
imported into the main EPIC interface, and the list of all materials associated with the spacecraft. 




| Thruster SPT100 Properties 


Name SPT100 
Type | Hall Thri 


x 0 


y 0.4 


z 0.4 


r 255 


g 192 


b 0 


dirX 0^ 
dirYpT 


Figure 1. Object ToolKit-generated spacecraft showing various entity definitions and the 

thruster (in this case a Stationary Plasma Thruster, SPT-100) location. 

Electric propulsion plumes 

Thruster plume definition in EPIC is performed using “PlumeTool” (Figure 3). The 2-D 
plume generation computational module allows the spacecraft designer to either import a plume 
map (produced by other codes) or generate it using the existing 2-D plume code. The 2-D plume 
code was written to model low-density EP exhausts such as those associated with electrostatic 
thrusters. The code uses parts from an existing, finite element (FE) code ( Gilbert ) that has been 
used in the past to solve a variety of problems involving Poisson’s equation (e.g. see Mikellides, 
I.G., et al , 7 ). The generated plume maps consist of the spatial distributions of ion density and 
velocity. The model consists of three main components: a Lagrangian algorithm for determining 
the expansion of the main ion beam, a Particle-in-Cell (PIC) solver for computing the dynamics 
of the charge-exchange plasma, and an ideal-scattering algorithm that calculates the density and 
velocity of ions that are produced by elastic collisions with neutral atoms. The physics and 
numerics of each component have been described in detail in previous publications. 1,8 ’ 9 ’ 1() 
Highlighted below are some of the fundamental principles behind the various plume model 
components. 







Figure 2. Main EPIC interface showing imported spacecraft object and thruster 
properties. 

2-D Plume Physics Models. The main ion beam is assumed to be a collisionless, singly- 
ionized, quasi-neutral plasma expanding under the influence of the electric field. The latter is 
assumed to exist due to the presence of density gradients in the plasma. By comparison to heavy- 
particle motion, electrons reach dynamic equilibrium at much smaller characteristic times. The 
electron inertia term may therefore be neglected in the equation of motion. In the absence of 
electron-ion collisions and magnetic fields, integration of the conservation of momentum for the 
electrons leads to the Boltzmann relation, which can be expressed in terms of the electric 
potential as follows: 

4>(n,T e ) = T e ln^/nj (1) 

where, <|) is the electric potential, n is the plasma density (n=ne=ni), T e is the electron 
temperature and a* is the reference plasma density at zero potential. 
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Hall-Effect Thruster 



Inputs Real Time Outputs Time-Averaged Outputs Details 
Simulation Time 4.000E 4 
Trial Potential 301.8 
Discharge Current 1.465 
Final Potential 300.3 
Exit Ion Current 1.334 
Exit Neutral Density 4.234E18 
Exit Plasma Density 1.940E17 
Exit Electron Temperature 9.235 
Exit Electric Field 1.516E4 
Exit Potential 194.6 
Exhaust Ion Current 1.449 
Exhaust Neutral Density 4.574E18 
Exhaust Plasma Density 7.812E16 
Average Ion Kinetic Energy 255.3 


Figure 3. PlumeTool’s Engine Specification window may be launched from the EPIC main 
interface. 

Ions are accelerated by the electric field, E = — V<|> according to, 

D Vj 

m i— = -e V< 1) (2) 

where m; and e are ion mass and electron charge, respectively. Since the drift velocity of the 
ions, Vj , is much greater than their thermal velocity, the high velocity ions are modeled as a 
fluid. The steady state, conservation equations of mass and momentum are solved in 2-D (R-Z) 
geometry. The numerical algorithm for computing the expansion of the main ion beam is based 
on a Lagrangian approach in which discrete mass elements, or otherwise “macroparticles” are 
released from the thruster exit and are then “tracked” using fundamental trajectory kinematics. 10 
The (normalized) current density profiles presently implemented along the exit of Ion and Hall 
thrusters are shown in Figure 4. For an Ion thruster the current density profile is assumed to vary 
parabolically with radius, i.e. as [l-(r/r out ET) ], where r 0U tET is the thruster outer radius. The ion 
velocity components V r and V z are determined by assuming that all ions are emitted from a point 




source behind the acceleration grid, with a maximum divergence angle of a (currently set at 20 
deg). Thus, 
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where Eb is the main beam energy. For Hall thrusters the current density and velocity profiles at 
the exit are based on extrapolation of integrated measurements from the Busek-Primex Hall 
thruster (BPT-4000). 5 



Figure 4. Normalized current density profiles for Ion (left) and Hall thrusters (right). 


The neutral gas density in space is mainly due to particles from the thruster and from the 
hollow cathode. The beam of neutrals from the thruster is computed using an annular anode gas 
flow model with isotropic emission from a ring. The profile of neutrals from the thruster is 
computed using two disk emissions defined by the solid angles subtended by each disk, and 
subtracting the smaller from the larger. 8 The flow rate of neutrals from the thruster exit is 
estimated using the known anode flow rate and propellant utilization. The hollow cathode is 
offset by a distance r H c from the thruster. Its axial location is assumed to be at z=0. The constant- 
temperature neutrals are emitted isotropically from the neutralizer, and their speed is estimated 
based on the same assumptions used for the thruster neutrals. 

Fast ions from the main beam undergo charge-exchange (CEX) with neutral particles, 
resulting in slow-moving ions and fast-moving neutrals, 


fast^ e + siowXe — > slow Xe + fast Xe 


( 4 ) 


In EPIC’s 2-D plume code, charge-exchange is computed using a two-dimensional PIC 
method. The rate of CEX-ion production rate h is determined by, 

n = n V:G n (5) 

CEX 1 CEX o 

In eqn (5) n is the particle density of neutrals and c is the cross section associated with 

O CEX 

the charge exchange collision. In contrast to the approach for the calculation of the main beam 
ions, where direct use of equation (1) provides the potential, the PIC algorithm solves the 2-D 
Poisson’s equation on a finite element grid and iterates until steady state CEX densities and 
potentials are consistent. The plasma density in Poisson’s equation is the sum of the main-beam 
and CEX densities. The first (computed by the Lagrangian method), and the prescribed neutral 
gas profile (above), are used as input for the calculation. The computed (total) ion particle 
density in Figure 5 shows clearly the CEX ion regions in the plume of an SPT-100. 
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Figure 5. Total ion particle density in the plume of an SPT-100 as computed by 
PlumeTool’s 2-D code. 

With the revival of nuclear sources in space flight many EP systems are projected to operate 
at high power levels (>10 kW) for increased specific impulse and thrust. Main beam ion energies 
in such missions can range in the thousands of electron-volts (>5000 eV) with ion flow rates 
exceeding a few milligrams per second (>3 mg/sec). In these propulsion systems the need to 
quantify the effects of high-energy ions (>300 eV) that are produced by elastic scattering 






between main beam ions and neutral particles becomes critical. The main reason is that these 
ions can be scattered to angles greater than those associated with the divergence of the main 
beam (which usually does not exceed 45 deg) and can therefore pose significant damage to 
surrounding spacecraft structures and diagnostics. An algorithm that determines the flux and 
energy of elastically scattered ions is part of PlumeTool; we outline here the basic approach. 

In the classical sense, the differential cross section associated with an atomic collision can be 
calculated from the deflection angle, 0(b,E c ,) n 


0(b,E c ) 


71 - 



dr 

-\/l — (b/r) 2 -(V(r)/E c ) 


( 6 ) 


where E c is the center-of-mass collision energy, b is the impact parameter, r is the interatomic 
distance, r m is the classical turning point (point of nearest approach), and V(r) is the interatomic 

potential energy. In the center of mass frame, G, is given by 0 = |®|, 0 < 0 < n . 


The differential cross section is obtained from equation (7) as follows: 
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( 7 ) 


where a is the collision cross section. The solid angle is given by dQ=27i sin0 d0. In the 
numerical model, given the angle and energy for each scattering pair (s-point to f-node as it will 
be shown later), the impact parameter is obtained by solving equation (6) using a Newton- 
Raphson method. The solution requires apriori knowledge of the interatomic potential function, 
V(r), which may in general have both attractive and repulsive contributions. At the high energies 
of interest in the present work (>300 eV) the deflection function barely exhibits a minimum as 
the scattering is almost solely governed by the repulsive part of the interaction potential. In such 
cases the classical approach (vs. the more rigorous quantum mechanical approach 12 ) is 
sufficiently accurate. Currently, only the repulsive potential function is included in the model. 
The coefficients that define the functional form of V(r) have been derived from averaged 
potentials by Amarouche, M., et al . 1 3 and are listed in Katz, I., et al. u The derivative in equation 
(7) is computed using a first order forward finite difference. 

We refer to Figure 6 to describe the formulation used for the determination of the density, 
average velocity and average energy of the elastically scattered ions. It is assumed that at the 
scattering point “s”, one scattering event occurs as the main beam ion flux, Fb=nb, s Ub, s impacts 
one stationary neutral particle. Once scattered, each ion does not undergo another collision. Then 
the particle flow rate dN/dt (#particles N per unit time t) through the incremental surface 
dS=p 2 dQ, is given by: 



where dQ is the solid angle subtended by dS at point “s”, Ij on (G) = \ 1(0) + j I(tt - 0) and p is the 
distance between “s” and field node “f The particle flow rate is also equal to, 

dN 

— = n sf u sf dS (9) 

dt 


where, n sf and u sf are the elastically scattered particle density and speed, respectively, at the 
field node “f Combining equations (8) and (9) we obtain, 


F _ n _4F b I IO „cos0 S f 

Kt ~ n sf u sf I 


(10) 


Equation (10) expresses the magnitude of the flux vector, F sf = |F sf | , at the field node location 
f (r f , z f ,0) due to scattering of the main beam flux Fb by a neutral at the scattering point 
s(r s , z s , (|)) . We note the additional term 4 cos 0 sl introduced in equation ( 1 0) by the 
transformation from the center-of-mass frame to the laboratory frame. 9 

ft is also assumed that the distance between nodes is much larger than the mean free path for 
ion-neutral collisions. Since both momentum and kinetic energy are conserved in the scattering 
event we can write, 

U sf 2 = U b,s 2 C °S 2 $sf (11) 


which also yields the density of scattered ions at “f ’, n sf as follows: 

_ F sf _ 4n b u b ,l, oil cos 9 sf _ 4n b I, 0l1 

11 sf ~~ ~~ 2 n _ 2 

u sf pu bjS cosO sf p 

The contributions to “f ’ from all points “s” along the scattering ring defined by each 
scattering node at s(r s ,z s ,0) , and by all nodes is summed as shown in equation (13), 

n f (r,z)=[^^n 0 dV 
v P 


( 12 ) 


(13) 


where no is the neutral particle density at the scattering node. 

Finally, we compute the average energy of the scattered ions at each field node using the 
energy flux and the particle flux as shown in equation (14) below. E S f is computed using equation 
( 11 ). 


E f (r,z) 


jE sf F sf n 0 dV 

V 

jF sf n 0 dV 

V 


(14) 




Figure 6. Elastic scattering geometry used to formulate the generalized algorithm in 

EPIC’s plume code. The four (red) corners define a portion of the plume computational 
plane. 

Comparisons with idealized test problems and measurements have been conducted for the 
purpose of validating the algorithm. The comparisons showed good agreement and are described 
in Gardner, B.M., et al., 10 and Mikellides, I.G., et al. 9 An example calculation for the BPT-4000 
plume is shown in Figure 7 (the computed main beam only is shown on the left). The computed 
particle density of the elastically scattered ions, when the thruster was operated in the laboratory, 
is depicted in Figure 7, right. The value of background pressure for this plume calculation was 
taken to be 3.3e-5 Torr. Figure 8 shows the plume map for an SPT-100 imported in the main 
EPIC interface. 



Figure 7. Computed main beam (left) and elastically-scattered ion density (right) (using 
the EPIC algorithm), in the plume of the BPT-4000 engine operating in the laboratory 
(P=3.3e-5 Torr). 



Figure 8. Imported SPT-100 plume showing ion particle density profile. 

1-D HET Code. In electric propulsion systems such as HETs for example, it is difficult to 
predict conditions at the exit using empirical measurements exclusively. Furthermore, the effect 
of thruster geometry and operating conditions on the exhaust, and consequently the surrounding 






S/C surfaces, is difficult to assess by simple extrapolation of measured quantities associated with 
a given thruster. Therefore, specifically for Hall thruster plumes EPIC incorporates a one- 
dimensional computer model of the acceleration region, which can be launched from within 
PlumeTool, as shown in Figure 3. The physics include collisions of electrons with neutrals, ions, 
and walls. The electron energy treatment includes ohmic heating, thermal diffusion and 
convection, and losses due to ionization, collisions with walls and atomic excitation. 15 

Plume-spacecraft interactions 

The ability to assess a variety of spacecraft interactions with electric propulsion plumes, in 
three dimensions, is the bulk of EPIC’s capability. The interactions models presently allowable 
by EPIC have been described in greater detail elsewhere. 16 We summarize here the most 
commonly used interactions capabilities of EPIC. It is noted apriori that EPIC does not account 
for any plume perturbations caused by the presence of the spacecraft. Therefore, except for 
straight-line shadowing, all interactions calculations are carried out using a fixed plume profile 
(as provided by PlumeTool). 

Fluxes to Surfaces. The ion flux f)jk (atoms m" 2 s' 1 ) at any point, i, on a surface j, due to 
plume component k, is calculated as follows: 

fijk = Pik v ik -n j (15) 

where p;k is the density at point “i” due to component “k”, v;k is the ion velocity (m/s) at point “i” 
of component “k”, and iij is the outward normal to surface “j”. Fluxes to points on surfaces 
account for the interference (“blocking”) by other spacecraft surfaces. Specifically, if a straight 
line between the point in question and the thruster orifice intercepts any other surface the flux is 
zero. 

2 1 

Surface Sputtering. The instantaneous sputtering rate (atoms m' s' ) of a spacecraft surface 
“j” at a point “i” due to thruster plume impingement, R S ij, is 

R | =Sk' Y ijk% 0«> 

where “k” is summed over the plume components and Yyt is the sputter yield at point “i” of 
surface “j” from plume component “k”. The sputter yield depends on the material, the energy of 
the ions impacting the surface and the angle between the flux vector and the surface-normal. All 
formulae are fits to sputtering measurements using the following functional form: 

Y(E, 0) = (a+bE)( 1.00-0.72 0+11.72 0 2 -3.13 0 3 -2.57 0 4 ) (17) 

1 7 

The yield is given in atoms/ion. The angular dependence is from Roussel et. al. 

Erosion/Deposition. The erosion/deposition rate at point “i” of surface “j” is the difference 
between the deposition and erosion rates. 



( 18 ) 
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If Rij > 0 it is a net deposition rate; if Ry < 0, it is a net erosion rate. 
The deposition rate is given by 


r?=i„ 


R»cos(e ljn )|fe 


(19) 


C 

where R n is average sputtering rate of all the points of surface “n,” Q m is the solid angle 

subtended by surface “n” viewed from point “i” and 0i jn is the angle between the normal of 
surface “j” and a ray from the centroid of surface “n” to the point “i”. 

The thickness or depth is calculated from the total particle deposition/erosion and the 
material “Density (gm cm" 3 )” and “Molecule mass” of the underlying material. Material 
properties are displayed on the “Material” sub-panel. The net erosion/deposition rate is computed 
at a randomly chosen selection of points over the mission (as specified in the “Monte Carlo” sub- 
panel) and then averaged. The mission averaged value is the product of the mission averaged rate 
and the “Mission Duration (days)” shown in the Orbit sub-panel. This ensures that time- 
dependent changes in spacecraft geometry (such as solar array sun-pointing) are included in the 
calculation. 

For example, one may be interested in determining how long the thruster can be operated 
before material from the solar array is eroded below a critical thickness. The user can plot 
erosion of a specific element as a function of time of thruster operation as shown in Figure 9. 

The solar array element chosen is shown on the spacecraft. Additional variables may be plotted 
in this manner. The green line depicted on the plot of Figure 9 is material deposited on the 
antenna element that is highlighted in red in Figure 2. 




Figure 9. Trade studies may be performed using the Plots/Trade Studies tab. The plot 
shows erosion and deposition on critical spacecraft surfaces as a function of thruster 
operation time. 

Induced Torques. Many electric thrusters produce torques that are comparable to or larger than 
those produced by gravity gradient, solar radiation, magnetic and aerodynamic effects, but are 
generally far smaller than torques resulting from chemical propulsion. Thrusters used for station- 
keeping or other purposes may produce unwanted torques, either directly or through the reflection of 
plume ions and neutrals from spacecraft surfaces. 

There are in general two limiting cases for the reflection of plume particles from surfaces. In 
specular elastic reflection the component of momentum normal to the surface is reversed. Reflection 
from rough surfaces leads to diffuse reflection in which the outgoing particle travels in a random 
direction, /.<?., particles are reflected in all directions regardless of the angle of incidence. In the case 
shown in Figure 10 (torque results are shown on the bottom left), specular reflection has been 
chosen. The general EPIC model used to determine the induced torques on the spacecraft during 
thruster operation accounts for contributions from the thrust and from the impingement of the 
exhaust on surfaces: 


r = SAR T x(-f T ) + SAR j xf j 
T j 


( 20 ) 


where, f j is the force imparted onto a surface j from plume particles, f T is the thrust vector, T is 
the torque, AR T is the position vector of a thruster from a reference point and ARj is the 
position vector of a surface j from a reference point. In the model, the contribution from plume 


impact is computed using one of the two limiting cases mentioned above: specular (elastic) 
reflection from the surface and fully diffuse reflection. 



Figure 10. Additional calculations and display options in EPIC. 

Summary 

The near and long-term benefits of electric propulsion are now widely recognized and 
accepted by both the government and the private sector. Solar-powered EP systems have been 
used successfully in numerous missions, with many more readied for near-term applications. 
Through Project Prometheus NASA now also focuses on nuclear electric propulsion (NEP) as 
the “way to defeat distance and time in exploring the outer planets.” However, integration of EP 
systems on spacecraft presents serious barriers to implementation. In the absence of accurate 
assessments, interactions between the propulsion system and the surrounding spacecraft 
components may threaten mission success. Quantifying such interactions within the design of 
specific spacecraft is a complex task with no comprehensive tools for the spacecraft engineer in 
existence. In response to the growing need by the electric propulsion community for a stand- 
alone design tool that addresses these issues, SAIC is developing the Electric Propulsion 
Interactions Code, EPIC. EPIC is a windows-based interactive, computer tool that allows the 
spacecraft designer to conduct the complete interactions process including geometry definition, 
EP plume generation, and interactions evaluation in 3-D. A fully functional beta version of the 
tool was delivered to the SEE Program on July, 2003. 
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